function [x,output]=agent12solverhog(x0,index,ewsb,waut);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table  ewdev statedev12 statedev13 statedev23
objval = [];
output=[];
options = optimset('LargeScale','off','Display','iter','TolX',0.0001);
[x,fval,exitflag]=fsolve(@objfun,x0,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);
d=index(3);
consu1 = x(1);
consu2 = x(2);
gammarfun=[x(3) x(4)];

wint=interpolation(scalefactor*ewsb, index,gammarfun);


output(1)=utility(rho,consu1)+beta*wint(1)/scalefactor;
output(2)=utility(rho,consu2)+beta*wint(2)/scalefactor;
output(3)=utility(rho,Y(j)-consu2-consu1)+beta*wint(3)/scalefactor;
f(1)=gammarfun(1) - consu2/consu1;
f(2)=gammarfun(2) - (Y(j)-consu2-consu1)/consu1;
f(3)= scalefactor*(utility(rho,S(j,1))+beta*ewdev(d,statedev12(j),1)'-output(1));
f(4)= scalefactor*(utility(rho,S(j,2))+beta*ewdev(d,statedev12(j),2)'-output(2));
end

end